{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "from matplotlib import animation\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import HTML\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, you will learn how to make use of the periodicity of the electrodes.\n",
    "\n",
    "As seen in [TB 4](../TB_04/run.ipynb) the transmission calculation takes a considerable amount of time. In this example, we will redo the *same* calculation, but speed it up (no approximations made).\n",
    "\n",
    "A large computational effort is made on calculating the self-energies which basically is inverting, multiplying and adding matrices, roughly 10–20 times per $k$-point, per energy point, per electrode.  \n",
    "For systems with large electrodes compared to the full device, this becomes more demanding than calculating the Green function for the system.  \n",
    "When there is periodicity in electrodes along the transverse semi-infinite direction (not along the transport direction) one can utilize Bloch's theorem to reduce the computational cost of calculating the self-energy.\n",
    "\n",
    "> In ***ANY*** calculation if you have periodicity, please ***USE*** it. \n",
    "\n",
    "In this example, you should scour the tbtrans manual on how to enable Bloch's \n",
    "theorem, and once enabled it should be roughly 3 - 4 times as fast, something that is non-negligible for large systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(orthogonal=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the below lines are differing from the same lines in [TB 4](../TB_04/run.ipynb), i.e. we save the electrode electronic structure *without* extending it 25 times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_elec = sisl.Hamiltonian(graphene)\n",
    "H_elec.construct(([0.1, 1.43], [0., -2.7]))\n",
    "H_elec.write('ELEC.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [TB 2](../TB_02/run.ipynb) for details on why we choose `repeat`/`tile` on the Hamiltonian object and not on the geometry, prior to construction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = H_elec.repeat(25, axis=0).tile(15, axis=1)\n",
    "\n",
    "center = H.geometry.center(what='xyz')\n",
    "\n",
    "H = H.remove(\n",
    "    H.geometry.close(center, R=10.)\n",
    ")\n",
    "\n",
    "dangling = [ia for ia in H.geometry.close(center, R=14.)\n",
    "                if len(H.edges(ia)) < 3]\n",
    "H = H.remove(dangling)\n",
    "edge = [ia for ia in H.geometry.close(center, R=14.)\n",
    "         if len(H.edges(ia)) < 4]\n",
    "edge = np.array(edge)\n",
    "\n",
    "# Pretty-print the list of atoms\n",
    "print(sisl.utils.list2str(edge + 1))\n",
    "\n",
    "H.geometry.write('device.xyz')\n",
    "H.write('DEVICE.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "Instead of analysing the same thing as in [TB 4](../TB_04/run.ipynb) you should perform the following actions to explore the available data-analysis capabilities of `sisl`. Please note the difference in run-time between example 04 and this example. Always use Bloch's theorem when applicable!\n",
    "\n",
    "*HINT* please copy as much as you like from example 04 to simplify the following tasks.\n",
    "\n",
    "1. Read-in the resulting file into a variable called `tbt`.\n",
    "2. In the following, we will concentrate on *only* looking at $\\Gamma$-point related quantities. I.e. all quantities should only be plotted for this $k$-point.  \n",
    "   To extract information for one or more $k$-points, you should look into the function:\n",
    "       \n",
    "       help(tbt.kindex)\n",
    "   \n",
    "   which may be used to find a resulting $k$-point index in the result file.\n",
    "   \n",
    "3. Plot the transmission ($\\Gamma$-point only). To extract a subset $k$-point you should read the documentation for the functions (*hint: `kavg` is the keyword you are looking for*).\n",
    "   - Full transmission\n",
    "   - Bulk transmission\n",
    "   \n",
    "4. Plot the DOS with normalization according to the number of atoms ($\\Gamma$ only)  \n",
    "   You may decide which atoms you examine.\n",
    "   - The Green function DOS\n",
    "   - The spectral DOS\n",
    "   - The bulk DOS\n",
    "   \n",
    "5. **TIME**: Do the same calculation using only tiling: `H_elec.tile(25, axis=0).tile(15, axis=1)` instead of `repeat`/`tile`. Which of `repeat` or `tile` are faster?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transmission"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Density of states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbt = sisl.get_sile('siesta.TBT.nc')\n",
    "# Easier manipulation of the geometry\n",
    "geom = tbt.geometry\n",
    "a_dev = tbt.a_dev # the indices where we have DOS\n",
    "# Extract the DOS, per orbital (hence sum=False)\n",
    "DOS = tbt.ADOS(0, sum=False)\n",
    "# Normalize DOS for plotting (maximum size == 400)\n",
    "# This array has *all* energy points and orbitals\n",
    "DOS /= DOS.max() / 400\n",
    "a_xyz = geom.xyz[a_dev, :2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "fig = plt.figure(figsize=(12,4));\n",
    "ax = plt.axes();\n",
    "scatter = ax.scatter(a_xyz[:, 0], a_xyz[:, 1], 1);\n",
    "ax.set_xlabel(r'$x$ [Ang]'); ax.set_ylabel(r'$y$ [Ang]');\n",
    "ax.axis('equal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If this animation does not work, then don't spend time on it!\n",
    "def animate(i):\n",
    "    ax.set_title('Energy {:.3f} eV'.format(tbt.E[i]));\n",
    "    scatter.set_sizes(DOS[i]);\n",
    "    return scatter,\n",
    "anim = animation.FuncAnimation(fig, animate, frames=len(tbt.E), interval=100, repeat=False)\n",
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learned lessons\n",
    "\n",
    "- Always use Bloch's theorem to greatly speed up the calculation speed\n",
    "- Extraction of number of coupled elements via `.edges` for the Hamiltonian.\n",
    "- Extending (i.e. `repeat`/`tile`) the Hamiltonian *after* it has been created (*very* fast!)\n",
    "- Extract data only for single $k$-points, the lesson learned is also applicable for a subset of $k$-points.\n",
    "- Extraction of various physical quantities from the `*.TBT.nc` file"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
